Propagation Delay Calculation and flight path distance estimation

In this tutorial we will see how to calculate the $t_0$ with the known sample. To ensure accurate estimation you should select your peak within the big mc range. You should also make your calculation for a small volume specially and for a small variation of voltage. Another pont is that you should choose a narrow window around your peaks.

# Activate intractive functionality of matplotlib
%matplotlib ipympl
# Activate auto reload 
%load_ext autoreload
%autoreload 2
%reload_ext autoreload
# import libraries
import os
import numpy as np
from sklearn import linear_model
import matplotlib.pyplot as plt
from ipywidgets import fixed, interact_manual, widgets
from ipywidgets import Output
from IPython.display import clear_output

# Local module and scripts
from pyccapt.calibration.calibration_tools import share_variables
from pyccapt.calibration.data_tools import data_tools, data_loadcrop, dataset_path_qt
from pyccapt.calibration.mc import mc_tools
from pyccapt.calibration.calibration_tools import widgets as wd
from pyccapt.calibration.calibration_tools import mc_plot

By clicking on the button below, you can select the dataset file you want to use. The dataset file can be in various formats, including HDF5, EPOS, POS, ATO, and CSV.

button = widgets.Button(
    description='load dataset',
)
@button.on_click
def open_file_on_click(b):
    """
    Event handler for button click event.
    Prompts the user to select a dataset file and stores the selected file path in the global variable dataset_path.
    """
    global dataset_path
    dataset_path = dataset_path_qt.gui_fname().decode('ASCII')
button

In case of recieving the error about pytable library, you have to install the pytables library with conda command. to do that you can open a new cell and copy the line below in it. Then just run it like other cells. The pytables library will be innstalled.

!conda install --yes --prefix {sys.prefix} pytables

From the dropdown lists below, you can select the instrument specifications of the dataset. The instrument specifications are the same as the ones used for the calibration process. Data mode is specify the dataset structure. The flight path length is the distance between the sample and the detector. The t0 is the time of flight of the ions with the lowest mass-to-charge ratio. The maximum mass-to-charge ratio is the maximum mass-to-charge ratio of tat you want to plot. You can also change it in te related cells. The detector diameter is the diameter of the detector.

tdc, pulse_mode, flightPathLength_d, t0_d, max_mc, det_diam = wd.dataset_instrument_specification_selection()
display(tdc, flightPathLength_d, max_mc)

# exctract needed data from Pandas data frame as an numpy array
# create an instance of the Variables opject
variables = share_variables.Variables()
variables.pulse_mode = pulse_mode

dataset_main_path = os.path.dirname(dataset_path)
dataset_main_path = os.path.dirname(dataset_main_path)
dataset_name_with_extention = os.path.basename(dataset_path)
variables.dataset_name = os.path.splitext(dataset_name_with_extention)[0]
variables.result_data_path = dataset_main_path + '/t_0_flight_path_distance/'
variables.result_data_name = variables.dataset_name
variables.result_path = dataset_main_path + '/t_0_flight_path_distance/'

if not os.path.isdir(variables.result_path):
    os.makedirs(variables.result_path, mode=0o777, exist_ok=True)
    
# Create data farame out of hdf5 file dataset
data = data_tools.load_data(dataset_path, tdc.value, mode='processed')
# extract data from the path and create the Variable object
data_tools.extract_data(data, variables, flightPathLength_d.value, max_mc.value)

The maximum time of flight: 5010

data

	x (nm)	y (nm)	z (nm)	mc_c (Da)	mc (Da)	high_voltage (V)	pulse	start_counter	t_c (ns)	t (ns)	x_det (cm)	y_det (cm)	pulse_pi	ion_pp
0	0.0	0.0	0.0	0.0	14.136714	5019.720215	1003.943970	3495	0.0	446.853577	2.964898	-0.169796	0	0
1	0.0	0.0	0.0	0.0	29.616535	5019.720215	1003.943970	3565	0.0	616.451904	-1.936327	0.088163	70	2
2	0.0	0.0	0.0	0.0	30.456534	5019.720215	1003.943970	4103	0.0	623.172729	-1.648980	0.672653	538	1
3	0.0	0.0	0.0	0.0	29.550233	5019.720215	1003.943970	4134	0.0	616.253052	-1.975510	-0.218776	31	1
4	0.0	0.0	0.0	0.0	28.966265	5019.720215	1003.943970	4205	0.0	605.431091	1.296327	-0.173061	71	1
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
10201104	0.0	0.0	0.0	0.0	141.314648	6347.270020	1269.453979	8768	0.0	1154.633423	-0.734694	-1.577143	403	1
10201105	0.0	0.0	0.0	0.0	29.461296	6347.270020	1269.453979	8799	0.0	548.825195	0.408163	1.508571	31	1
10201106	0.0	0.0	0.0	0.0	29.600126	6347.270020	1269.453979	9443	0.0	547.803345	-1.165714	0.097959	485	1
10201107	0.0	0.0	0.0	0.0	29.719453	6347.270020	1269.453979	9771	0.0	555.038513	-1.645714	1.296327	328	1
10201108	0.0	0.0	0.0	0.0	29.823040	6347.270020	1269.453979	10327	0.0	556.574707	-0.982857	1.933061	556	1

10201109 rows × 14 columns

Extract the data from the dataset file in numpy array format.

# exctract needed data from Pandas data frame as an numpy array

# exctract needed data from Pandas data frame as an numpy array
pulse_mode.value = 'voltage'
variables.dld_high_voltage = data['high_voltage (V)'].to_numpy()
variables.dld_pulse = data['pulse'].to_numpy()
variables.dld_t = data['t (ns)'].to_numpy()
variables.dld_x_det = data['x_det (cm)'].to_numpy()
variables.dld_y_det = data['y_det (cm)'].to_numpy()
# calculate the mc based on t_0 = 0
variables.mc = mc_tools.tof2mc(variables.dld_t, 0, variables.dld_high_voltage, variables.dld_x_det, variables.dld_y_det, flightPathLength_d.value, 
                                         variables.dld_pulse, mode=pulse_mode.value)

interact_manual(data_loadcrop.plot_crop_experiment_history, data=fixed(data), variables=fixed(variables), max_tof=widgets.FloatText(value=variables.max_tof), frac=widgets.FloatText(value=1.0),
                bins=fixed((1200,800)), figure_size=fixed((7,3)),
               draw_rect=fixed(False), data_crop=fixed(False), pulse_plot=widgets.Dropdown(options=[('False', False), ('True', True)]), dc_plot=widgets.Dropdown(options=[('True', True), ('False', False)]),
                pulse_mode=widgets.Dropdown(options=[('voltage', 'voltage'), ('laser', 'laser')]), save=widgets.Dropdown(options=[('True', True), ('False', False)]),
               figname=widgets.Text(value='exp_hist'));

C:\Users\APTUser\AppData\Local\Temp\ipykernel_10944\2178511811.py:1: DeprecationWarning: on_submit is deprecated. Instead, set the .continuous_update attribute to False and observe the value changing with: mywidget.observe(callback, 'value').
  interact_manual(data_loadcrop.plot_crop_experiment_history, data=fixed(data), variables=fixed(variables), max_tof=widgets.FloatText(value=variables.max_tof), frac=widgets.FloatText(value=1.0),

plot the histogram of the mass spectrum.

interact_manual(
        mc_plot.hist_plot,
        variables=fixed(variables),
        bin_size=widgets.FloatText(value=0.1),
        log=widgets.Dropdown(options=[('True', True), ('False', False)]),
        target=widgets.Dropdown(options=[('mc', 'mc'), ('tof', 'tof')]),
        mode=widgets.Dropdown(options=[('normal', 'normal'), ('normalized', 'normalized')]),
        prominence=widgets.IntText(value=100),
        distance=widgets.IntText(value=100),
        percent=widgets.IntText(value=50),
        selector=fixed('peak'),
        figname=widgets.Text(value='hist'),
        lim=widgets.IntText(value=variables.max_mc),
        peaks_find_plot=widgets.Dropdown(options=[('True', True), ('False', False)]),
        peaks_find=fixed(True),
        range_plot=fixed(False),
        plot_ranged_ions=fixed(False),
        ranging_mode=fixed(False),
        selected_area_specially=fixed(False),
        selected_area_temporally=fixed(False),
        save_fig=widgets.Dropdown(options=[('True', True), ('False', False)]),
        print_info=fixed(True),
        figure_size=fixed((9, 5)));

C:\Users\APTUser\AppData\Local\Temp\ipykernel_10944\966379802.py:1: DeprecationWarning: on_submit is deprecated. Instead, set the .continuous_update attribute to False and observe the value changing with: mywidget.observe(callback, 'value').
  interact_manual(

For the selected peak you have to select the mass-to-charge ratio of the peak. You can select the mass-to-charge ratio from the dropdown list below. You can also select the charge of the peak. After that you can add the peak to the list by clicking on the add button. You can also delete the selected peak from the list by clicking on the delete button. You can reset the list by clicking on the reset button.

isotopeTableFile = '../../../files/isotopeTable.h5'
dataframe = data_tools.read_hdf5_through_pandas(isotopeTableFile)
elementsList = dataframe['element']
elementIsotopeList = dataframe['isotope']
elementMassList =  dataframe['weight']
abundanceList = dataframe['abundance']

elements = list(zip(elementsList, elementIsotopeList, elementMassList, abundanceList))
dropdownList = []
for element in elements:
    tupleElement = ("{} ({}) ({:.2f})".format(element[0], element[1], element[3]), "{}({})[{}]".format(element[0], element[1], element[2]))
    dropdownList.append(tupleElement)

chargeList = [(1,1,),(2,2,),(3,3,),(4,4,)]
dropdown = wd.dropdownWidget(dropdownList,"Elements")
dropdown.observe(wd.on_change)


chargeDropdown = wd.dropdownWidget(chargeList,"Charge")
chargeDropdown.observe(wd.on_change_charge)

wd.compute_element_isotope_values_according_to_selected_charge()

buttonAdd = wd.buttonWidget("ADD")
buttonDelete = wd.buttonWidget("DELETE")
buttonReset = wd.buttonWidget("RESET")

display(dropdown, chargeDropdown, buttonAdd, buttonDelete, buttonReset)

def buttonAdd_f(b,):
    with out:
        clear_output(True)
        wd.onClickAdd(b, variables)
        display()
def buttonDelete_f(b,):
    with out:
        clear_output(True)
        wd.onClickDelete(b, variables)
        display()
def buttonResett_f(b,):
    with out:
        clear_output(True)
        wd.onClickReset(b, variables)
        display()

listMaterial = buttonAdd.on_click(buttonAdd_f)
buttonDelete.on_click(buttonDelete_f)
buttonReset.on_click(buttonResett_f)
# listMaterial = buttonAdd.on_click(wd.onClickAdd)
# buttonDelete.on_click(wd.onClickDelete)
# buttonReset.on_click(wd.onClickReset)
out = Output()
display(out)

# The correct peak location based on your sample
peaks_ideal = variables.list_material
print('peak ideal mass:', peaks_ideal)
peak_x = variables.peaks_x_selected
print('peak actual mass:', peak_x)
print('The peak index are:', variables.peaks_index_list)

peak ideal mass: [1.01, 13.49, 26.98]
peak actual mass: [1.0024705770121665, 14.435576305489002, 29.272140841418643]
The peak index are: [0, 1, 2]

Here for each peak a mask will be created.

pick_ions_plot = []
mask = np.zeros(len(variables.mc), dtype=bool)
mc_seb_ideal = np.zeros(len(variables.mc))
# creat mask for each peak base on the peak loc. and window size
index = 0
for i in variables.peaks_index_list:
    mask_tmp = np.logical_and((variables.x_hist[round(variables.peak_widths[2][i])] < variables.mc), (variables.mc < variables.x_hist[round(variables.peak_widths[3][i])]))
    # Count the number of True values in the original mask_tmp
    true_indices = np.where(mask_tmp)[0]
    num_true_values = len(true_indices)

    # If there are more than 5000 True values, randomly choose 5000 of them
    if num_true_values > 5000:
        selected_indices = np.random.choice(true_indices, size=5000, replace=False)

        # Create a new mask with the same length as the original, initialized with False
        new_mask = np.full_like(mask_tmp, False)

        # Set the selected indices to True in the new mask
        new_mask[selected_indices] = True

        # Update the original mask_tmp with the new mask
        mask_tmp = new_mask
    
    print('peak left and right sides are:', variables.x_hist[round(variables.peak_widths[2][i])], variables.x_hist[round(variables.peak_widths[3][i])])
    mask = np.logical_or(mask, mask_tmp)
    mc_seb_ideal[mask_tmp==True] = peaks_ideal[index]
    index += 1
    bb = np.zeros(len(variables.mc))
    
    pick_ions_plot.append(mask_tmp)

peak left and right sides are: 1.0024705770121665 1.1027176346873666
peak left and right sides are: 14.134835132463401 14.836564536189803
peak left and right sides are: 28.570411437692243 29.873623187469846

Here we plot ions in each peak base on the TOF and (x,y) position

for i in range(len(pick_ions_plot)):
    fig1, ax1 = plt.subplots(figsize=(6, 6))
    dld_x_masked = variables.dld_x_det[pick_ions_plot[i]]
    dld_y_masked = variables.dld_y_det[pick_ions_plot[i]]
    dld_t_masked = variables.dld_t[pick_ions_plot[i]]
    if len(dld_x_masked) > 1000:
        mask_plot = np.random.randint(0, len(dld_x_masked), 1000)
        x = plt.scatter(dld_x_masked[mask_plot], dld_t_masked[mask_plot], color="blue", label='X', alpha=0.1)
        y = plt.scatter(dld_y_masked[mask_plot], dld_t_masked[mask_plot], color="red", label='Y', alpha=0.1)
    else:
        x = plt.scatter(dld_x_masked, dld_t_masked, color="blue", label='X', alpha=0.1)
        y = plt.scatter(dld_y_masked, dld_t_masked, color="red", label='Y', alpha=0.1)
    ax1.set_ylabel("Time of flight (ns)", color="red", fontsize=10)
    ax1.set_xlabel("position (cm)", color="red", fontsize=10)
    plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
    plt.legend(handles=[x, y], loc='upper right')
    plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
    plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
    plt.savefig(variables.result_path + 'position_peak.png', format="png", dpi=600)
    plt.savefig(variables.result_path + 'position_peak.svg', format="svg", dpi=600)
    plt.show()

As you saw the TOF changes base on the (x,y) of the events. Therefore we creat a mask to only select the ions in center (8m m*8mm) of detector. This helps to cansel out the variation in TOF due to hit position.

Reformulate the equation:

$$t = d(\sqrt{\frac{\frac{m}{n}}{2eV}})-t_{0} $$

# use mask_equal to have equal number of ions for each peak
# only peak the value in the center of detector 10mm * 10mm
detector_squre = 0.5
vol_variation = 5500 # peak only for 200 voltage variation for vol_variation +/- 100 

fig1, ax1 = plt.subplots(figsize=(6, 6))
dld_x_masked = variables.dld_x_det[mask]
dld_y_masked = variables.dld_y_det[mask]
dld_t_masked = variables.dld_t[mask]
dld_highVoltage_masked = variables.dld_high_voltage[mask]
dld_pulse_masked = variables.dld_pulse[mask]

mask_tmp_middle = np.logical_and((np.abs(dld_x_masked) < detector_squre), (np.abs(dld_y_masked) < detector_squre))
mask_tmp_mvoltage = mask_tmp_mvoltage = np.logical_and((dld_highVoltage_masked < (np.mean(dld_highVoltage_masked)+200)), (dld_highVoltage_masked > (np.mean(dld_highVoltage_masked)-200)))
mask_tmp_mvoltage = np.ones(len(mask_tmp_middle), dtype=bool)
mask_f = np.logical_and(mask_tmp_mvoltage, mask_tmp_middle)

dld_x_masked = dld_x_masked[mask_f]
dld_y_masked = dld_y_masked[mask_f]
dld_t_masked = dld_t_masked[mask_f]
dld_highVoltage_masked = dld_highVoltage_masked[mask_f]
dld_pulse_masked = dld_pulse_masked[mask_f]

mc_seb_reg_masked = mc_seb_ideal[mask]
mc_seb_reg_masked = mc_seb_reg_masked[mask_f]

x = plt.scatter(dld_x_masked, dld_t_masked, color="blue", label='X', alpha=0.1)
y = plt.scatter(dld_y_masked, dld_t_masked, color="red", label='Y', alpha=0.1)
ax1.set_ylabel("Time of flight (ns)", color="red", fontsize=10)
ax1.set_xlabel("position (cm)", color="red", fontsize=10)
plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
plt.legend(handles=[x, y], loc='upper right')
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
plt.savefig(variables.result_path + 'center.png' , format="png", dpi=600)
plt.savefig(variables.result_path + 'center.svg' , format="svg", dpi=600)

plt.show()

We calculate the ${t_0}$ base on:

$$t_{0} = \frac{\sum_{n=1}^n{\left (t - L_{flight} \sqrt{\frac{m/n}{2eV}} \right )}}{n}$$

seb_t = dld_t_masked * 1E-9  # tof in s
                                                               
seb_factor = np.sqrt(mc_seb_reg_masked * 1.66E-27 / (2 * 1.6E-19 * dld_highVoltage_masked))

seb_factor = seb_factor * 1E6
seb_t = seb_t * 1E9

t0_seb_fixed = np.mean(np.array([seb_t]).squeeze(0) - (flightPathLength_d.value * np.array([seb_factor]).squeeze(0).reshape(-1, 1)))
print('Linear fixed path lenght -- the flight path lenght(slop): {:.2f}'.format(flightPathLength_d.value), '(mm)', '\nthe corrected t_0(intercept): {:.2f}'.format(t0_seb_fixed), '(ns)')

Linear fixed path lenght -- the flight path lenght(slop): 110.00 (mm) 
the corrected t_0(intercept): -13.83 (ns)

# Plot outputs
# fig1, ax1 = plt.subplots(figsize=(5.5/2.54, 5.5/2.54))
fig1, ax1 = plt.subplots(figsize=(5.5, 5.5))
peaks_data = plt.scatter(seb_factor, seb_t, label='Ions', alpha=1, color='slategray')
axes = plt.gca()
x_vals = np.array(axes.get_xlim())

linear_fix, = plt.plot(x_vals, t0_seb_fixed + flightPathLength_d.value * x_vals, '--', color='red', label='line' )

plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
# plt.legend(handles=[linear_fix, peaks_data], loc='lower right')

ax1.set_ylabel("Time of Flight (ns)", fontsize=8)

ax1.set_xlabel(r"$\sqrt{\frac{\frac{m}{n}}{2e \alpha V}} (ns/mm)$", fontsize=8)
plt.tight_layout()
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
plt.savefig(variables.result_path + 'fixed_flight_path.svg', format="svg", dpi=600)
plt.savefig(variables.result_path + 'fixed_flight_path.png', format="png", dpi=600)

plt.show()

linear = linear_model.LinearRegression()
linear.fit(np.array([seb_factor]).squeeze(0).reshape(-1, 1), np.array([seb_t]).squeeze(0))
d_seb = linear.coef_.item()
t0_seb = linear.intercept_.item()

print('Linear -- the corrected flight path lenght(slop): {:.2f}'.format(d_seb), '(mm)', '\nthe corrected t_0(intercept): {:.2f}'.format(t0_seb), '(ns)')

Linear -- the corrected flight path lenght(slop): 106.98 (mm) 
the corrected t_0(intercept): -2.58 (ns)

# Plot outputs

# fig1, ax1 = plt.subplots(figsize=(5.5/2.54, 5.5/2.54))
fig1, ax1 = plt.subplots(figsize=(5.5, 5.5))
peaks_data = plt.scatter(seb_factor, seb_t, color="slategray", label='Ion', alpha=0.1)
axes = plt.gca()
x_vals = np.array(axes.get_xlim())


linear, = plt.plot(x_vals, t0_seb + d_seb * x_vals, '--', color='r', label='fit' )


plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)

# plt.legend(handles=[peaks_data, linear, rigid, huber,lasso], loc='lower right')
ax1.set_ylabel("Time of flight (ns)", fontsize=8)
ax1.set_xlabel(r"$\sqrt{\frac{\frac{m}{n}}{2e \alpha V}} (ns/mm)$", fontsize=8)
plt.legend(handles=[linear, peaks_data], loc='lower right')
plt.tight_layout()
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])

plt.savefig(variables.result_path + 'regression.svg', format="svg", dpi=600)
plt.savefig(variables.result_path + 'regression.png', format="png", dpi=600)
    

plt.show()

Plot the m/c with new ${t_0}$:

# recalulate the mc based on the new t_0
mc = mc_tools.tof2mc(variables.dld_t, 14, variables.dld_high_voltage, variables.dld_x_det, variables.dld_y_det, flightPathLength_d.value, 
                                         variables.dld_pulse, mode=pulse_mode.value)

mc_hist = mc_plot.AptHistPlotter(mc[mc < 40], variables)
y, x = mc_hist.plot_histogram(bin_width=0.1, label='mc', steps='stepfilled', log=True, grid=True, fig_size=(9, 5))
peaks, properties, peak_widths, prominences = mc_hist.find_peaks_and_widths(prominence=100, distance=10, percent=50)
mc_hist.plot_peaks()
mc_hist.plot_hist_info_legend(label='mc', bin=0.1, background=None, loc='right')